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Abstract —The most efficient signal edge-preserving smoothing 
filters, e.g., for denoising, are non-linear. Thus, their acceleration 
is challenging and is often performed in practice by tuning 
filter parameters, such as by increasing the width of the local 
smoothing neighborhood, resulting in more aggressive smoothing 
of a single sweep at the cost of increased edge blurring. We 
propose an alternative technology, accelerating the original filters 
without tuning, by running them through a special conjugate 
gradient method, not affecting their quality. The filter non¬ 
linearity is dealt with by careful freezing and restarting. Our 
initial numerical experiments on toy one-dimensional signals 
demonstrate 20x acceleration of the classical bilateral filter and 
3-5x acceleration of the recently developed guided filter. 

Index Terms —conjugate gradient algorithm, edge-preserving 
denoising, low-pass filters 

L Introduction 

This papeiQ is concerned with noise removal from a given 
noisy signal, which is a basic problem in signal processing, 
with many applications, e.g., in image denoising Q. Modern 
denoising algorithms preserve signal details while removing 
most of the noise. A very popular denoising filter is the 
bilateral filter (BE), which smooths signals while preserving 
edges, by taking the weighted average of the nearby pixels. 
The weights depend on both the spatial distance between the 
sampling locations and similarity between signal values, thus 
providing local adaptivity to the input signal. Bilateral filtering 
has initially been proposed in ifT^ as an intuitive tool without 
theoretical justification. Since then, connections between BE 
and other well-known filtering techniques such as anisotropic 
diffusion, weighted least squares, Bayesian methods, kernel 
regression and non-local means have been explored; see, e.g., 
survey El. 

We make use of the graph-based framework for signal 
analysis developed in El, Q, where polynomial low-pass 
filters based on the BE coefficients are proposed. A nice 
introduction to signal processing on graphs is found in m. 

A single application of BE can be interpreted as a vertex 
domain transform on a graph with pixels as vertices, intensity 
values of each node as the graph signal, and filter coefficients 
as link weights that capture the similarity between nodes. The 
BE transform is a special nonlinear anisotropic diffusion, cf. 
El, El, determined by the entries of the graph Laplacian 
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matrix, which are related to the BE weights. The eigenvectors 
and eigenvalues of the graph Laplacian matrix allow us to 
extend the Eourier analysis to the graph signals or images as 
in H and perform frequency selective filtering operations on 
graphs, similar to those in traditional signal processing. 

Another very interesting smoothing filter is the guided filter 
(GE), recently proposed in El, 0, and included into the 
MATLAB image processing toolbox. Some ideas behind GE 
are developed in lH. According to our limited experience, GE 
is faster than BE. The authors of El advocate that GE is 
gradient preserving and avoids gradient reversal artifacts in 
contrast to BE, which is not gradient preserving. 

The smoothing explicit filters similar to BE and GE can be 
interpreted as matrix power iterations, which are, in general 
case, nonlinear, or equivalently, as explicit integration in time 
of the corresponding nonlinear anisotropic diffusion equation 
El, El. The suitable graph Laplacian matrices are determined 
by means of the graph-based interpretation of these power 
iterations. Our main contribution is accelerating the smoothing 
filters by means of a special variant of the conjugate gradient 
(CG) method, applied to the corresponding graph Laplacian 
matrices. To avoid oversmoothing, only few iterations of the 
CG acceleration can be performed. We note that there exist 
several nonlinear variants of the CG algorithm, see, e.g.. 
However, the developed theory is not directly applicable in our 
case because it is not clear how to interpret the vector L{x)x 
as a gradient of a scalar function of the signal x, where L{x) 
is a graph Laplacian matrix depending on a signal x. 

IT Bilateral eilter (BE) 

We consider discrete signals defined on an undirected graph 
Q = (V, 8), where the vertices V = {1, 2,..., W} denote, e.g., 
time instances of a discrete-time signal or pixels of an image. 
The set of edges £ = {{i,j)} contains only those pairs of 
vertices i and j that are neighbors in some predefined sense. 
We suppose in addition that a spatial position pi is assigned 
to each vertex i S V so that a distance \\pi —pj\\ is determined 
between vertices i and j. 

Let a;[y], j G V, be a discrete function, which is an input 
signal to the bilateral filter. The output signal y[i] is the 
weighted average of the signal values in x[f\: 
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The weights Wij are defined for {i,j) £ £ in terms of a 
guidance signal g[i]: 


Wij = exp 
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where Od and ar are the filter parameters 02 . The guidance 
signal g\i\ is chosen depending on the purpose of filtering. 
When g coincides with the input x, the bilateral filter is 
nonlinear and called self-guided. 

The weights Wij determine the adjacency matrix W of 
the graph G. The matrix W is symmetric, has nonnega¬ 
tive elements and diagonal elements equal to 1. Let D be 
the diagonal matrix with the nonnegative diagonal entries 
di = ■ Thus, the BF operation ([U is the vector 

transform defined by the aid of the matrices W{g) and D{g) 
as y = D~^Wx = x — D~^Lx, where L = D — W \s called 
the Laplacian matrix of the weighted graph G with the BF 
weights. The eigenvalues of the matrix D~"^W are real. The 
eigenvalues corresponding to the highest oscillations lie near 
the origin. 

The BF transform y = D~^Wx can be applied iteratively, 
(i) by changing the weights Wij at each iteration using the 
result of the previous iteration as a guidance signal g, or (ii) 
by using the fixed weights, calculated from the initial signal 
as a guidance signal, for all iterations. The former alternative 
results in a nonlinear filter. The latter produces a linear filter, 
which may be faster, since the BF weights are computed only 
once in the very beginning. 

An iterative application of the BF matrix transform is the 
power iteration with the amplification matrix D~^W. Slow 
convergence of the power iteration can be boosted by the aid 
of suitable Krylov subspace iterative methods Q, HI. 


III. Guided filter (GF) 


Algorithm 1 Guided Filter (GF) 

Input: X, g, p, e 

Output: y 

TnedTlg = fmeaniy^ p') 
meaUx = fmean{x,p) 

COTTg = fmeanig- * g,p) 

COrVgx = fmeanig- * P) 
vavg = corvg — meaug. * meaug 
covgx = corvgx — meaUg. * meaUx 
a = coVgx-Hvarg -f e) 
b = meanx — a. * meaUg 

TTLCdTla — fmeani^-> P) 

IJlCdTlfj — fmeanib^ P^ 

y = mearia- * g + medUb 

Algorithm 1 is a pseudo-code of GF proposed in HU, where 
X and y are, respectively, the input and output signals on the 
graph G, described in section |II] GF is built by means of a 
guidance signal g, which equals x in the self-guided case. 
The function fmeani-,p) denotes a mean filter of a spatial 
radius p. The constant e determines the smoothness degree 
of the filter—the larger e the larger smoothing effect. The 


dot preceded operations .* and ./ denote the componentwise 
multiplication and division. A typical arithmetical complexity 
of the GF algorithm is 0{N), where N is the number of 
elements in x, see El. 

The guided filter operation of Algorithm 1 is the matrix 
transform y — Wig)x, where the implicitly constructed 
transform matrix W{g) has the following entries, see El : 
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The mean filter fmeani', p) is applied in the neighborhoods Wk 
of a spatial radius p around all vertices fc £ V. The number of 
pixels in Wk is denoted by |tLi|, the same for all k. The values 
Pk and are the mean and variance of g over ojk- The matrix 
W is symmetric and satisfies the property = 1. 

The standard construction of the graph Laplacian matrix 
gives L = I — W, because di — J^j — 1’ i.^. 
matrix D is the identity. The eigenvalues of L(p) are real 
nonnegative with the low frequencies accumulated near 0 and 
high frequencies near 1. Application of a single transform 
y = Wx attenuates the high frequency modes of x while 
approximately preserving the low frequency modes, cf. Q, 
0 . Similar to the BF filter, the guided filter can be applied 
iteratively. When the guidance signal g is fixed, the iterated 
GF filter is linear. When g varies, for example, g = x for the 
self-guided case, the iterated GF filter is nonlinear. 


IV. Conjugate gradient acceleration 

Since the graph Laplacian matrix L is symmetric and 
nonnegative definite, the iterative application of the transform 
y = D~^WX can be accelerated by adopting by the CG 
technology. We use two variants of CG: 1) with the fixed 
guidance equal to the input signal or to the clean signal, 2 ) 
with the varying guidance equal to the current value of x. 

Algorithm 2 Truncated PCG(fcniax) 

Input: xo, g, k^ax Output: x 
X = Xo; r = Wig)x — Dig)x 
for k — \ ^^ kxaax 1 do 
s = D~^ig)r; 7 = s'^r 

if fc = 1 then p = s else /3 = p = s + fip 

endif 

q = D{g)p - Wig)p; a = 7 /(p^q) 

X = X + ap; r = r — aq; "fold = 7 

endfor 

Algorithm 2 is the standard preconditioned conjugate gra¬ 
dient algorithm formally applied to the system of linear 
equations Lx = 0 and truncated after fcmax evaluations of the 
matrix-vector operation Lx. The initial vector xq is a noisy 
input signal. This variant of the CG algorithm has first been 
suggested in 0 . 

Algorithm 3 is a special nonlinear preconditioned CG with 
^max restarts, formally applied to L{x)x = 0 and truncated 
after fcmax iterations between restarts. Restarts are necessary 
because of nonlinearity of the self-guided filtering. 









Algorithm 3 Truncated PCG(fcniax) with /max restarts 


PSNR = 20.0494, SNR = 13.0753 


Input- Xq, /t^max? ^max Output- X 
X = Xo 

for / = 1, . . . , /max do 
r = W{x)x — D{x)x 
for k = 1, , fcmax - 1 do 

s = D~^{x)r', 7 = s'^r 

if /c = 1 then p = s else /3 = "fl^oid', p = s + f3p 

endif 

q = D{x)p — W (x)p; a = {p'^q) 

X = X -\- ap; r = r — aq; 'fold = 7 

endfor 

endfor 


V. Numerical experiments 

As a proof of concept, our MATLAB tests use the clean 
1-dimensional signal Xc of length N = 4730 shown in Fig¬ 
ure [T] We choose this rather difficult, although 1 -dimensional, 
example to better visually illustrate both the denoising and 
edge-preserving features of the filters. The noisy signal, also 
displayed in Figure [T] is the same for all tests and given by 
the formula xq = Xc + rj, where a Gaussian white noise rj 
has zero mean and variance cr^ = 0.01. The bilateral filter is 
used with ad = 0.5 and = 0.1. The neighborhood width 
in BF equals 5 so that the band of W consists of 5 diagonals. 
The guided filter is used with e = 0.001 and the neighborhood 
width 3. The matrix W of GF also has 5 diagonals. 

The CG accelerated BF/GF is called CG-BF/CG-GF. Typi¬ 
cal numerical results of the average performance are displayed. 

The signal error after denoising is x — xq, where x stands 
for the output denoised signal. We calculate the peak signal-to- 
noise ratio (PSNR) and signal-to-noise ration (SNR). The pa¬ 
rameters are manually optimized to reach the best possible 
match of the signal errors in the compared filters, resulting in 
indistinguishable error curves in our figures. 

The results in Figures |2] and |3] are obtained by the iterated 
BF and GF filters with the fixed guidance g = Xc and by 
CG-BF and CG-GF implemented in Algorithm 2 with the 
same fixed guidance g = Xc- These tests are performed only 
for comparison reasons because the clean signal guidance Xc 
seems to be ideal for the best possible denoising results. 

We say that Algorithm 3 uses /max x fcmax iterations, if 
it executes /max restarts with the fcmax evaluations L(x)x 
between restarts. The best denoising performance for our test 
problem is achieved with the following iteration combinations 
of the self-guided CG-BF; 31 x 3, 17 x 4, 12 x 5, 9 x 6, 7 x 7, 
6 X 8, 5 X 9, 4 X 10, 3 X 11, 2 X 19. The best combinations 
for the self-guided CG-GF are 11 x 3, 7 x 4, 5 x 5, 4 x 6, 
3x7. Figures |4] and |5] show the results after 3 x 11 iterations 
of CG-BF and 5x5 iterations of CG-GF. 

The numerical tests demonstrate about 20-times reduction 
of iterations for the self-guided bilateral filter and 3-times 
reduction of iterations for the guided filter with self-guidance 
after the conjugate gradient acceleration. It is also interesting 
to observe that both filters with the properly chosen parameters 
and iteration numbers produce almost identical output signals. 



Fig. 1. Clean and noisy signals. 
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Fig. 2. 500 BF iterations versus 20 CG-BF iterations with the guidance Xc- 
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Fig. 3. 90 GF iterations versus 13 CG-GF iterations with the guidance Xc- 


VI. Conclusion 

Iterative application of BF and GF, including their non¬ 
linear self-guided variants, can be drastically accelerated by 
using CG technology. Our future work concerns developing 
automated procedures for choosing the optimal number of CG 
iterations and investigating CG acceleration for 2D signals. 
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Fig. 4. 600 iterations of the self-guided BF versus 3 x 11 iterations of CG-BF. 
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Fig. 5. 75 iterations of the self-guided GF versus 5x5 iterations of CG-GF. 
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